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We exactly solve a Fokker-Planck equation by determining its eigenvalues and eigenfunctions: we 
construct nonlinear second-order differential operators which act as raising and lowering operators, 
generating ladder spectra for the odd and even parity states. These are staggered: the odd-even 
separation differs from even-odd. The Fokker-Planck equation describes, in the limit of weak damp- 
ing, a generalised Ornstein-Uhlenbeck process where the random force depends upon position as 
well as time. Our exact solution exhibits anomalous diffusion at short times and a stationary 
non-Maxwellian momentum distribution. 

PACS numbers: 05.40.-a,05.45.-a,02.50.-r,05.40.Fb 



There are only few physically significant systems with 
ladder spectra (exactly evenly spaced energy levels). Ex- 
amples are the harmonic oscillator and the Zeeman- 
splitting Hamiltonian. In this letter we introduce a pair 
of exactly solvable eigenvalue problems. We use them 
to solve two (closely related) physically significant ex- 
tensions of a classic problem in the theory of diffusion, 
namely the Ornstein-Uhlenbeck process. These eigen- 
value problems have ladder spectra. Our systems differ 
from the examples above in that their spectra consist 
of two ladders which are staggered; the eigenvalues for 
eigenfunctions of odd and even symmetry do not inter- 
leave with equal spacings. We introduce a new type of 
raising and lowering operators in our solution, which are 
nonlinear second-order differential operators. Our gen- 
eralised Ornstein-Uhlenbeck systems exhibit anomalous 
diffusion at short times, and non-Maxwellian velocity dis- 
tributions at equilibrium; we obtain exact expressions 
which are analogous to results obtained for the standard 
Ornstein-Uhlenbeck process. 

1. Ornstein- Uhlenbeck processes. Before we discuss our 
extension of the Ornstein-Uhlenbeck process, we describe 
its usual form. This considers a particle of momentum 
p subjected to a rapidly fluctuating random force f(t) 
and subject to a drag force — jp, so that the equation of 
motion is p — —~/p + /(£)• The random force has statis- 
tics (/(*)) = 0, {f(t)f(t')) = C(t - t') (angular brackets 
denote ensemble averages throughout). If the correlation 
time r of f(i) is sufficiently short (jt <C 1), the equation 
of motion may be approximated by a Langevin equation: 
dp = — 7pdt + dw, where the Brownian increment dw 
has statistics (dw) — and (dw 2 ) — 2D$dt. The diffu- 
sion constant is D = \ di(/(i)/(0)). This problem 
is discussed in many textbooks (for example Q|); it is 
easily shown that the variance of the momentum (with 
the particle starting at rest) is 

(p 2 (t)) = [l-exp(-27t)]A>/7, (1) 

that the equilibrium momentum distribution is Gaussian 
and that the particle (of mass m) diffuses in space with 
diffusion constant T> x — Do/m 2 -/ 2 . 



In many situations the force on the particle will be 
a function of its position as well as of time. Our ex- 
tension of the Ornstein-Uhlenbeck process is concerned 
with what happens in this situation when the damping 
is weak. We consider a force f(x,t) which has mean 
value zero, and a correlation function (f(x,t)f(x',t')) — 
C(x — x' , t — t'). The spatial and temporal correlation 
scales are £ and r respectively. If the momentum of the 
particle is large compared to pq = m^/r, then the force 
experienced by the particle decorrelates more rapidly 
than the force experienced by a stationary particle. Thus, 
if the damping 7 is sufficiently weak that the particle is 
accelerated to a momentum large compared to po, the 
diffusion constant characterising fluctuations of momen- 
tum will be smaller than Do- The impulse of the force 
on a particle which is initially at x = in the time from 
t = to t = At is 

Aw = dt f(pt/m,t) + 0(At 2 ) . (2) 
Jo 

If At is large compared to r but small compared to I/7, 
we can estimate (Aw 2 ) = 2D(p)At, where 

D{p) = \f dtC(pt/m,t). (3) 

When p ^ po we recover D(p) = Dq. When p 3> po, we 
can approximate to obtain 

D(p) = -^l +0 (p- 2 ) ,D X = — dX C(X,0) . (4) 

\P\ 2p J-oo 

If the force is the gradient of a potential, f(x,t) — 
dV(x,t)/dx, then D\ = 0. In this case, expanding the 
correlation function in Q in its second argument gives 
D(p) ~ D 3 Pq/\p 3 \, where D 3 may be expressed as an 
integral over the correlation function of V(x, t). To sum- 
marize: the momentum diffusion constant is a decreas- 
ing function of momentum, such that D(p) ~ H^ 1 for 
a generic random force, or D(p) ~ H~ 3 for a gradient 
force. 
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2. Fokker- Planck equation. The probability density for 
the momentum, P(p,t), satisfies a Fokker-Planck (gener- 
alised diffusion) equation. Following arguments in pj, 
one obtains 



d t P = d p ( 1P P + D{ P )d p P) 



(5) 



Sturrock Q introduced a related Fokker-Planck equation 
(without the damping term) and also gave an expression 
for D{p) analogous to equation J2J. His work was applied 
to the stochastic acceleration of particles in plasmas, and 
subsequent contributions have concentrated on refining 
models for the calculation of D(p) (see, for example, 0, 
4]). In the following we obtain exact solutions to (J5J 
in the cases where D(p) = Dipo/\p\ (which we consider 
first) and D(p) = D^p^/\p\ 3 (which can be treated in the 
same way, and for which we discuss the results at the end 
of the paper). 

Introducing dimensionless variables (t' — jt and z — 
Pil/PoDi) 1 / 3 ), the Fokker-Planck equation for the case 
where D(p) ocl/p becomes 



d t >P = dJzP 



1 d z P) = FP. 



(6) 



It is convenient to transform the Fokker-Planck operator 
F to a Hermitian form, which we shall refer to as the 
Hamiltonian operator 



H = P, 



l_\zl,d_J_d_ () 
2 4 dz \z\ dz { 1 



where Pq(z) oc exp(— |z| 3 /3) is the stationary solution 
satisfying FPq = 0. We solve the diffusion problem by 
constructing the eigenfunctions of the Hamiltonian op- 
erator. In the following we make free use of the Dirac 
notation pj of quantum mechanics to write the equations 
in a compact form and to emphasise their structure. 

3. Summary of principal results. We start by list- 
ing our results (for the case of random forcing, where 
D(p) oc l/|p|). We construct the eigenvalues A„ and 
eigenfunctions ip n {z) of the operator H. We identify rais- 
ing and lowering operators A + and A which map one 
eigenfunction to another with respectively one more or 
one fewer nodes. We use these to show that the spectrum 
of H consists of two superposed equally spaced spectra 
(ladder spectra) for even and odd parity states: 



A+ = -3t 



A,; = -(3n + 2) , n = 0. 



(8) 



The spectrum of the Hamiltonian J7J) is displayed on the 
rhs of of Fig. 1. It unusual because the odd-even step is 
different from the even-odd step, due to the singularity 
of the Hamiltonian at z = 0. Our raising and lowering 
operators allow us to obtain matrix elements required for 
calculating expectation values. One interesting result is 
the variance of the momentum for a particle starting at 
rest at t = 0: 



(P 2 (t)) = 



p 0J DiA 2/3 3 7 / 6 r(2/3) 



2tt 



(1 



_ e -37«)2/3 



(9) 
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FIG. 1: The spectrum of H (right) is the sum of two equally 
spaced (ladder) spectra and \„ shifted w.r.t each other 
(left). 



This is reminiscent of equation Q for the stan- 
dard Ornstein-Uhlenbeck process, however exhibits 
anomalous diffusion for small times. At large times 
(p 2 {t)) converges to the expectation of p 2 with the sta- 
tionary (non-Maxwellian) momentum distribution 



P (p) =Mexp(- 1 \p\ 3 /(p D 1 )) 



(10) 



(TV is a normalisation constant). At large times the dy- 
namics of the spatial displacement is diffusive (x 2 (t)) ~ 
2T> r t with diffusion constant 



(paD^ 3 7T3- 5 / 6 (I 1 2 5 5 
m 2 7 5/3 2F(2/3) 2 32 V3' 3' 3' 3' 3' 



(11) 



(here -F32 is a hypergeometric function). At small times, 
by contrast, we obtain anomalous diffusion 

(x 2 (*)) = C x (bo#i) 2/ V 8/3 ™- 2 ) t 8/3 (12) 

where the constant C x is given by ^ below. 

4. Ladder operators and eigenfunctions. The eigen- 
function of the Fokker-Planck equation JJJJ are alter- 
nately even and odd functions, defined on the inter- 
val (—00,00). The operator H, describing the limit- 
ing case of this Fokker-Planck operator, is singular at 
z = 0. We identify two eigenfunctions of H by in- 
spection, ipo(z) = Cq exp(— |z| 3 /6) which has eigenvalue 

value A 



and ipQ (z) = C z\z\ exp(— |z| 3 /6) with eigen- 
-2. These eigenfunctions are of even and 







odd parity, respectively (zero and one node, respectively). 
Our approach to determining the full spectrum will be to 
define a raising operator A + which maps any eigenfunc- 
tion (z) to its successor with the same parity, Vvi-iO 2 ) > 
having two additional nodes. 

We now list definitions of the operators we use: these 
include the definitions of the raising and lowering opera- 
tors, A + and A, as well as an alternative representation 
of the Hamiltonian: 



a ± = {d z ±z\z\/2) 



A -- 
H 



a \z\ 
- a~\z\ 



A+ 
G 



a \z\ 



(13) 
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FIG. 2: Shows (p 2 (t)} and (x 2 (t)}. Computer simulation 
of the equations of motion p — —"{p + f(x, i) and mx = p 
(symbols); theory, eqs. 1101 and 126H . red lines. Also shown 
are the limiting behaviours for (x 2 (t)}, lilt and 1121 . at long 
and short times (dashed lines). In the simulations, C(X, t) = 
cr 2 exp[-X 2 /(2£ 2 ) - t 2 /(2T 2 )}. The parameters were m = 1, 
7 = 1CT 3 , £ = 0.1, r = 0.1, and a = 20. 



Note that A + is the Hcrmitian conjugate of A. The most 
significant properties of the operators A and A + are 



[H, A] = 3i and [H, A+] = SA^ 



(14) 



(the square brackets are commutators). These expres- 
sions show that the action of A and A + on any eigenfunc- 
tion is to produce another eigenfunction with eigenvalue 
increased or decreased by three, or else to produce a func- 
tion which is identically zero. The operator A + adds two 
nodes, and repeated action of A + on i/)q(z) and i^q{z) 
therefore exhausts the set of eigenfunctions. Together 
with Xq = and = — 2 this establishes that the spec- 
trum of H is indeed iJSJ. Some other useful properties of 
the operators of equation (|13fl are 



A] = 3(H+G) , H-G = I , A+A = H 2 +2H . (15) 

We represent the eigenfunctions by of H by kets \tp~) 
and \ipn)- The actions of A and A + are 

A + \^t) = Ct +1 \^t + i) and A\i>t) = C^t-i) ■ 

(16) 

where (using JTSJ) we have C± = y/in(in^ 2). 

A peculiar feature of A and A + is that they are of sec- 
ond order in d/dz, whereas other examples of raising and 
lowering operators are of first order in the derivative. The 
difference is associated with the fact that the spectrum 
is a staggered ladder: only states of the same parity have 
equal spacing, so that the raising and lowering operators 
must preserve the odd-even parity. This suggests replac- 
ing a first order operator which increases the quantum 
number (total number of nodes) by one with a second 
order operator which increases the quantum number by 
two, preserving parity. 

There is an alternative approach to generating the 
eigenfunctions of (J5J ■ This equation falls into one of the 
classes considered in ;6j , and we have written down first- 
order operators which map one eigenfunction into an- 
other. However, these operators are themselves functions 



of the quantum number n, making the algebra cumber- 
some. We have not succeeded in reproducing our results 
with the 'Schrodinger factorisation' method. 

5. Propagator and correlation functions. The propa- 
gator of the Fokker-Planck equation dfK — FK can be 
expressed in terms of the eigenvalues A^ and cigcnfunc- 
tions <W=P " 1/2 C(z) otF: 

OO 

K(y, z; = E E e MKt') ■ (17) 

n=0 <r=±l 

Here y is the initial value and z is the final value of 
the coordinate. The expansion coefficients a^(y) are de- 
termined by the initial condition K(y, z;0) = S(z — y), 

— 1/2 

namely a°(y) = P fpniu)- l n terms of the eigenfunc- 
tions of H we have 

K(y,z;t') = Y,Pv 1 '\v)€Ay)P 1 J\zWn(z)eMKt') ■ 

Tier 

(18) 

The propagator determines correlation functions. As- 
suming zq = we obtain for the expectation value of 
a function O(z) at time t 

(O(z(<0)>-E7®(^ + |O(z)|^)exp(A+t') . (19) 

Similarly for the correlation function of 0(z(t2)) and 

o(z(h)) (with t' 2 > t[ > o) 

{0(z(t' 2 ))0(z(t[))) = J2 ff7^(^0 + |O(i)IC) 

x «|0(z)|V+) ex P K(t' 2 - t[) + A+i'r] . (20) 

6. Momentum diffusion. To determine the time- 
dependence of (p 2 {t)) we need to evaluate the ma- 
trix elements Yo n — (-0^~|z 2 |?/>+). A recursion for 
these elements is obtained as follows. Let Fon+i = 
(^\z 2 A+\iJ+)/C+ +1 . Write zA+ = zG + z(A+ - G) = 
z(H + z{A+ - G). It follows 

(^\z 2 A + \ri) = (\+-i)Y 0n +(^\z 2 (A+-GM+) • (21) 

Using (A + — G) = —za~ and [l 3 ,a~] = — 3z 2 we ob- 
tain Fo n+ i = A+ + 21on/Cn+i) and together with Y 0Q = 
3 7 / 6 r(2/3)/(27r) this gives 



V ' x/2W 2 (3n-2) y/Tfc + 1) V ' 
We also find 

and after performing the sum in l|19fl we return to dimen- 
sional variables. The final result is ©. 
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7. Spatial diffusion. The time-dependence of (x 2 ) is 
determined in a similar fashion, using l(20|) and ((23(1 : 

= J_ 3 _L /"dti /d^ . (24) 

The matrix elements Z mn = are found by a 

recursion method, analogous to that yielding J22J): 



thus evaluate the small-t' behaviour by approximating 
the sums in ((26(1 by integrals and using the asymptotic 
form for the coefficients Am- There is a non-integrable 
divergence of A(k,l) as k — > I, which can be cancelled 
by using the sum rule in equation 128(1 . We obtain the 
limiting behaviour ((12(1 with 



-C 



z „ 



( _ ir -n^! (m + n + 1)r(2/3) (25) 
D7T 



y/Y{n + l)r(m + 1/3) T(n - m + 1/3) 



dy 



dx x~ s ' 3 
a(x) — a(xy) 



(30) 



1 



y 



xa'(x) 



y2/3(l_y)5/3 



y/T(m + l)r(n + 5/3) T(n - m + 2) 
for Z > m — 1 and zero otherwise. Using (|20|l 
(po^i) 2/3 



TO 2 7 5 / 3 



with A fei - (V^(0)/^ o + (0))Z o; Z fei and 



where a(x) — [1 — exp(— x)]/a; and C = 
3~ 7/3 r(2/3)/(27r 2 ). The integral is convergent and 
can be evaluated numerically to give C x — 0.57 . . .. This 
is in good agreement with a numerical evaluation of the 
sum 1(26(1. as shown in Fig. [21 
53 53 ^ T fci(0 ( 26 ) 8. Gradient-force case. In the case where the force 

is the gradient of a potential function, we have D(p) = 
D 3 Pq/\p\ 3 + 0(p~ 4 ) 7}. In dimensionless variables the 
Fokker-Planck equation is 



oo oo 



T kl (t') = [ dti / dt' 2 e x i 
Jo Jt{ 

+ f dt[ [ tl dt' 2 e^^- t ^ +x t^. (27) 
Jo Jo 

We remark upon an exact sum-rule for the A)~i , and also 
on their asymptotic form for k 1, I 1: 



&,P = d,(zP+ \z\- 3 d z P) = FP 



(31) 



5> w = , A,-g 2/3)2 



k + l 



k=0 



3 l/3 47r 2 fc2/3;4/3(; _ fc )5/3 



(28) 



We now show how to derive the limiting behaviours 111(1 
and I12|) . shown as dashed lines in Fig. [21 At large 
time x evolves diffusively: (x 2 ) ~ 2T> x t, with the diffu- 
sion constant obtained from the equilibrium momentum 
correlation function: 



instead of ©. This Fokker-Planck equation has 
the non-Maxwellian equilibrium distribution Pq(z) — 
exp(— |z| 5 /5). The raising and lowering operators are 
of the form A + = a~\z\~ 3 a~~ and A = a + \z\~ 3 a + 
with a ± = (<9 Z ± z\z\ 3 /2). The analogue of (JTUl is 
= 5A, [H,A + ] = -5A+, and the eigenvalues 
are 0, —4, —5, —9, —10, —14,-15, .... In this case, too, a 
closed expression for example for (p 2 (t)} can be obtained 
which is analogous to l|TJl. but which exhibits anomalous 
diffusion: 



(P 2 t) = 



p a D 3 \ 2/5 5 2 / 5 sin(7r/5) 



-57*^2/5 



(32) 



x r(3/5)r(4/5)(l 



2m 2 7 V 7 / ^-^oo 



/oo 
dt' (zTZt'+r) 
-OO 



-1 /poA 



2to 2 7 \ 7 



2/3 



7 2 



A,, 



(29) 



which evaluates to 1(11(1 . At small values of i' the dou- 
ble sum 1(26(1 is dominated by the large-fc, I terms. We 



The short-time anomalous diffusion is consistent with the 
asymptotic scaling obtained in 0, Q for the special case 
of undamped stochastic acceleration. 

A full account of our results, and their extension to 
higher dimensions, will be published elsewhere. 
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